当clusterProfiler遇见stringdb...
上次搞的ProjectYulab，还是有小伙伴响应的。当前最好的一个案例就是我说支持PPI网络的，直接上截图。
趁着出来徐州开会，晚上醒醒酒的时间，把这些代码整理了。
那么在当前github版本的clusterProfiler中，我们把这块功能给加进去了。直接上实例吧。
require(aplot)
require(ggplot2)
require(clusterProfiler)
genes <- c("TP53", "BRCA1", "BRCA2", "S100A6")
g1 <- getPPI(genes, taxID="9606")
plot(g1)
一条plot指令就出图了。当然这只是igraph默认的画图，我们留待以后再加强一下画图的功能。
新加入的getPPI函数，你的输入只需要基因向量，以及相应的物种ID即可。
而这个物种ID，我们也搞了一个函数getTaxID来帮助你获取。
> getTaxID("Homo sapiens")
[1] "9606"
那么你肯定要问，你想把一些相互作用的蛋白也纳入进去呢？这有何难，用参数add_nodes指令数量就可以了。
g2 <- getPPI(genes, taxID='9606', add_nodes=10)
plot(g2)
接下来，你肯定要问，那和富集分析有什么关系？
我们来跑个例子吧。
data(geneList, package="DOSE")
de = names(geneList)[1:100]
x = enrichKEGG(de)
y = setReadable(x, 'org.Hs.eg.db', 'ENTREZID')
然后我们拿第一条富集的结果出来看。
> y[1,]
ID Description GeneRatio BgRatio pvalue p.adjust
hsa04110 hsa04110 Cell cycle 12/58 157/8398 5.077653e-10 6.14396e-08
qvalue
hsa04110 5.719041e-08
geneID
hsa04110 CDC45/CDC20/CCNB2/NDC80/CCNA2/CDK1/MAD2L1/CDT1/TTK/AURKB/CHEK1/TRIP13
Count
hsa04110 12
假如我就想看看富集于某个通路的基因其编码的蛋白质之间是个什么样的相互作用，那么现在就很简单了。富集分析的结果直接可以做为getPPI的输入，然后可以传一下富集分析中的第几条通路，是想要看的。
g = getPPI(y, 1)
plot(g)
要加蛋白进入，同样指定一下数目，就加进去了。
g2 = getPPI(y, 1, add_nodes=10)
plot(g2)
那么你还想问，能不能几条通路的基因一起进去做PPI网络，这当然可以。
聪明的你可能已经发现，为什么这些PPI网络都这么高度连通？因为stringdb就是这样，它默认是会把功能相关的（比如同表达）也做给连上边，如果我们就只要相互作用，那么边会少一些的，这个可能通过参数设定。
g4 = getPPI(y, 1:3, network_type = "physical")
plot(g4)
最后的最后，我们再配合一个aplot，那么我们就可以拼出来不同的通路相应的PPI了。
pp <- list(
~plot(getPPI(y, 1, network_type = "physical")),
~plot(getPPI(y, 2, network_type = "physical"))
)
names(pp) <- y$Description[1:2]
pg <- plot_list(gglist = pp, ncol=1)
pg
当然需要强调一下的是，我们还可以自己搞点基于ggplot2画网络图的函数，然后画图会更好一些。
好了，再回到开头，然后可以有更多的小伙伴加入到ProjectYulab中来，练练手还是不错的。
写完了，感觉还不错，酒也醒了，喝酒和写bug更配。写完bug该洗洗睡了，明早赶飞机。